

### This file produces Figures I1 through I4 that reproduce the main analysis
### for Hypothesis 3 based on the observed value approach. It requires the 
### output from the script "Main analyses"

library(plyr)
library(boot)
library(MASS)
library(tidyverse)
library(DAMisc)

##### Figure I1 on Hypothesis 1 #####

# function to predict probabilities
avg.pred.fun.gen <- function(model, invlink, varnames, values, sims, probs){
  obs.data <- model$data
  samp.coef <- t(mvrnorm(sims,coef(model),vcov(model)))
  hyp.data <- do.call(expand.grid,values)
  avg.preds <- adply(hyp.data,1,function(x){
    obs.data[,varnames] <- x
    des.mat <- model.matrix(model$formula,data = obs.data)
    all.preds <- invlink(des.mat %*% samp.coef)
    quantile(colMeans(all.preds),probs)
  })
  avg.preds
}

# set general parameter
invlink <- inv.logit

### for extraversion

# specify parameters
varnames <- c("extraST")
values <- list(extraST = c(-1, 1))

# predict
predictH1a_extra_oc <- data.frame(avg.pred.fun.gen(mH1_set1_2, invlink = invlink, varnames = varnames, 
                 values = values, sims = 250, probs = c(0.025, 0.5, 0.975)))
colnames(predictH1a_extra_oc) <- c("x", "conf.low", "predicted", "conf.high")

# plot
pH1a_extra_oc <- ggplot(predictH1a_extra_oc) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("-1 sd", "+1 sd")) +
  xlab("extraversion") +
  ylab("predicted probability of activism") +
  ylim(0.35, 0.65) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH1a_extra_oc

### for gender

# specify parameters
varnames <- c("female")
values <- list(female = c(0, 1))

# predict
predictH1a_female_oc <- data.frame(avg.pred.fun.gen(mH1_set1_2, invlink = invlink, varnames = varnames, 
                                       values = values, sims = 250, probs = c(0.025, 0.5, 0.975)))
colnames(predictH1a_female_oc) <- c("x", "conf.low", "predicted", "conf.high")

# plot
pH1a_female_oc <- ggplot(predictH1a_female_oc) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("male", "female")) +
  xlab("gender") +
  ylab("") +
  ylim(0.35, 0.65) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH1a_female_oc


### for Moscow

# specify parameters
varnames <- c("moscow")
values <- list(moscow = c(0, 1))

# predict
predictH1a_moscow_oc <- data.frame(avg.pred.fun.gen(mH1_set1_2, invlink = invlink, varnames = varnames, 
                                       values = values, sims = 250, probs = c(0.025, 0.5, 0.975)))
colnames(predictH1a_moscow_oc) <- c("x", "conf.low", "predicted", "conf.high")

# plot
pH1a_moscow_oc <- ggplot(predictH1a_moscow_oc) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("no", "yes")) +
  xlab("lives in Moscow") +
  ylab("") +
  ylim(0.35, 0.65) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH1a_moscow_oc

ggsave("Figure_I1.pdf", gridExtra::arrangeGrob(pH1a_extra_oc, pH1a_female_oc, pH1a_moscow_oc, ncol = 3),
       width = 7, height = 3)


##### Figure I2 on Hypothesis 2 #####

predictH2b_extra_ov <- data.frame(ordAveEffPlot(mH2_set1_4, "extraST", data = data2, plot = F))

# plot
predictH2b_extra_ov$outcome <- c(rep("none", 25), rep("sporadic", 25), rep("regular", 25))
colnames(predictH2b_extra_ov) <- c("predicted", "conf.low", "conf.high", "cat", "extraST", "outcome") 

predictH2b_extra_ov$outcome <- factor(predictH2b_extra_ov$outcome, 
                                      levels = c("none", "sporadic", "regular"))

p.H2_polr_ov <- ggplot(predictH2b_extra_ov) +
  geom_smooth(aes(x=extraST, y=predicted),
              size=.5, se = F, color = "black") +
  geom_ribbon(aes(x=extraST, ymin = conf.low, ymax = conf.high), 
              alpha = .2) +
  facet_wrap(factor(predictH2b_extra_ov$outcome)) +
  xlab("extraversion") +
  ylab("predicted probability of campaign work intensity") +
  theme_bw(); p.H2_polr_ov

ggsave("Figure_I2.pdf", width = 9, height = 4.5)


##### Figure I3 (first part) for Hypothesis 3 (within 2013 survey) #####

### for extraversion

# specify parameters
varnames <- c("extraST")
values <- list(extraST = c(-1, 1))

# predict
predictH3a_extra_oc <- data.frame(avg.pred.fun.gen(mH3_set1_2, invlink = invlink, varnames = varnames, 
                                                   values = values, sims = 250, probs = c(0.025, 0.5, 0.975)))
colnames(predictH3a_extra_oc) <- c("x", "conf.low", "predicted", "conf.high")

# plot
pH3a_extra_oc <- ggplot(predictH3a_extra_oc) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("-1 sd", "+1 sd")) +
  xlab("extraversion") +
  ylab("predicted probability of being\nin opp. activist group") +
  ylim(0.4, 0.8) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3a_extra_oc

### for agreeableness

# specify parameters
varnames <- c("agreeST")
values <- list(agreeST = c(-1, 1))

# predict
predictH3a_agree_oc <- data.frame(avg.pred.fun.gen(mH3_set1_2, invlink = invlink, varnames = varnames, 
                                                   values = values, sims = 250, probs = c(0.025, 0.5, 0.975)))
colnames(predictH3a_agree_oc) <- c("x", "conf.low", "predicted", "conf.high")

# plot
pH3a_agree_oc <- ggplot(predictH3a_agree_oc) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("-1 sd", "+1 sd")) +
  xlab("agreeableness") +
  ylab("\n") +
  ylim(0.4, 0.8) +
  #ggtitle("") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3a_agree_oc


### for gender

# specify parameters
varnames <- c("female")
values <- list(female = c(0, 1))

# predict
predictH3a_female_oc <- data.frame(avg.pred.fun.gen(mH3_set1_2, invlink = invlink, varnames = varnames, 
                                                   values = values, sims = 250, probs = c(0.025, 0.5, 0.975)))
colnames(predictH3a_female_oc) <- c("x", "conf.low", "predicted", "conf.high")

# plot
pH3a_female_oc <- ggplot(predictH3a_female_oc) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("male", "female")) +
  xlab("gender") +
  ylab("\n") +
  ylim(0.4, 0.8) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3a_female_oc


### for Moscow

# specify parameters
varnames <- c("moscow")
values <- list(moscow = c(0, 1))

# predict
predictH3a_moscow_oc <- data.frame(avg.pred.fun.gen(mH3_set1_2, invlink = invlink, varnames = varnames, 
                                                    values = values, sims = 250, probs = c(0.025, 0.5, 0.975)))
colnames(predictH3a_moscow_oc) <- c("x", "conf.low", "predicted", "conf.high")

# plot
pH3a_moscow_oc <- ggplot(predictH3a_moscow_oc) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("no", "yes")) +
  xlab("lives in Moscow") +
  ylab("\n") +
  ylim(0.4, 0.8) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3a_moscow_oc

ggsave("Figure_I3a.pdf", gridExtra::grid.arrange(pH3a_extra_oc, pH3a_agree_oc, pH3a_female_oc, pH3a_moscow_oc, 
                               ncol = 4, top = "2013 survey"),
       width = 10, height = 3)


##### Figure I3 (second part) for Hypothesis 3 (with pooled survey) #####

### for extraversion

# specify parameters
varnames <- c("extraST")
values <- list(extraST = c(-1, 1))

# predict
predictH3b_extra_oc <- data.frame(avg.pred.fun.gen(mH3_set1_4, invlink = invlink, varnames = varnames, 
                                                   values = values, sims = 250, probs = c(0.025, 0.5, 0.975)))
colnames(predictH3b_extra_oc) <- c("x", "conf.low", "predicted", "conf.high")

# plot
pH3b_extra_oc <- ggplot(predictH3b_extra_oc) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("-1 sd", "+1 sd")) +
  xlab("extraversion") +
  ylab("predicted probability of being\nin opp. activist group") +
  ylim(0.15, 0.65) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3b_extra_oc

### for agreeableness

# specify parameters
varnames <- c("agreeST")
values <- list(agree = c(-1, 1))

# predict
predictH3b_agree_oc <- data.frame(avg.pred.fun.gen(mH3_set1_4, invlink = invlink, varnames = varnames, 
                                                   values = values, sims = 250, probs = c(0.025, 0.5, 0.975)))
colnames(predictH3b_agree_oc) <- c("x", "conf.low", "predicted", "conf.high")

# plot
pH3b_agree_oc <- ggplot(predictH3b_agree_oc) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("-1 sd", "+1 sd")) +
  xlab("agreeableness") +
  ylab("\n") +
  ylim(0.15, 0.65) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3b_agree_oc


### for gender

# specify parameters
varnames <- c("female")
values <- list(female = c(0, 1))

# predict
predictH3b_female_oc <- data.frame(avg.pred.fun.gen(mH3_set1_4, invlink = invlink, varnames = varnames, 
                                                    values = values, sims = 250, probs = c(0.025, 0.5, 0.975)))
colnames(predictH3b_female_oc) <- c("x", "conf.low", "predicted", "conf.high")

# plot
pH3b_female_oc <- ggplot(predictH3b_female_oc) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("male", "female")) +
  xlab("gender") +
  ylab("\n") +
  ylim(0.15, 0.65) +
  #ggtitle("") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3b_female_oc


### for Moscow

# specify parameters
varnames <- c("moscow")
values <- list(moscow = c(0, 1))

# predict
predictH3b_moscow_oc <- data.frame(avg.pred.fun.gen(mH3_set1_4, invlink = invlink, varnames = varnames, 
                                                    values = values, sims = 250, probs = c(0.025, 0.5, 0.975)))
colnames(predictH3b_moscow_oc) <- c("x", "conf.low", "predicted", "conf.high")

# plot
pH3b_moscow_oc <- ggplot(predictH3b_moscow_oc) +
  geom_pointrange(aes(x = factor(x), y = predicted,
                      ymin = conf.low, ymax = conf.high)) +
  geom_text(aes(x = factor(x), y = predicted, label = round(predicted, 2)),
            size=4, hjust = 1.3) +
  scale_x_discrete(labels = c("no", "yes")) +
  xlab("lives in Moscow") +
  ylab("\n") +
  ylim(0.15, 0.65) +
  #ggtitle("") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)); pH3b_moscow_oc

ggsave("Figure_I3b.pdf", gridExtra::grid.arrange(pH3b_extra_oc, pH3b_agree_oc, pH3b_female_oc, pH3b_moscow_oc, 
                                  ncol = 4, top = "pooled 2013/Navalny"),
       width = 10, height = 3)


##### Figure I4 for interaction #####

# specify parameters
varnames <- c("extraST", "agreeST")
values <- list(extraST = seq(min(data1$extraST), max(data1$extraST),
                             length.out = 50),
               agreeST = c(-1,1))

# predict for 2013 survey
predictH4a_oc <- data.frame(avg.pred.fun.gen(m4_set1_1, invlink = invlink, varnames = varnames, 
                                                    values = values, sims = 250, probs = c(0.025, 0.5, 0.975)))
colnames(predictH4a_oc) <- c("x", "group", "conf.low", "predicted", "conf.high")


# predict for pooled sample
predictH4b_oc <- data.frame(avg.pred.fun.gen(m4_set1_2, invlink = invlink, varnames = varnames, 
                                             values = values, sims = 250, probs = c(0.025, 0.5, 0.975)))
colnames(predictH4b_oc) <- c("x", "group", "conf.low", "predicted", "conf.high")


# join both
predictH4_oc <- rbind(predictH4a_oc, predictH4b_oc)
predictH4_oc$sample <- c(rep("2013 survey", 100), rep("pooled: 2013/Navalny", 100))

# plot
predictH4_oc$group <- factor(predictH4_oc$group, levels = c("-1", "1"))

ggplot(predictH4_oc) +
  geom_smooth(aes(x=x, y=predicted, color = factor(group)),
              size=.5, se = F) +
  geom_ribbon(aes(x=x, ymin = conf.low, ymax = conf.high, fill = factor(group)), 
              alpha = .2) +
  scale_color_discrete(name="agreeableness",
                       labels = c("-1 sd", "+1 sd")) +
  scale_fill_discrete(name="agreeableness",
                      labels = c("-1 sd", "+1 sd")) +
  xlab("extraversion") +
  ylab("predicted probability of being\nin opp. activist group") +
  facet_wrap(~sample, 
             ncol = 3) +
  theme_bw()

ggsave("Figure_I4.pdf", width = 7, height = 4) 




